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PREFACE 


The work described in this report was performed by the Infor- 
mation Systems Division of the Jet Propulsion Laboratory. 
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THE COVARIANCE MATRIX FOR THE SOLUTION VECTOR 
OF AN EQUALITY -CONSTRAINED LEAST 
SQUARES PROBLEM 

1. INTRODUCTION 

Consider the linear least squares problem 

Ex Sf 

subject to the linear equality constraints 

Cx = d 

We refer to this as Problem LSE denoting Lease Squares with Eqxiality 
constraints. 

Methods for solving Problem LSE are described in Chapters 20 - 22 
of Ref. (1). 

In this note we describe methods for computing the covariance matrix 
V for the solution vector x. Different methods of computing V will be discussed 
which are convenient for use with each of the different solution algorithms 
given in Ref. (1), Any reference to a Chapter, Section, or Page without further 
qualification is to be understood to refer to Ref. (1). 

We assume throughout that the covariance matrix of f is the identity 
matrix. If the covariance matrix of f is known to be something other than the 
identity matrix then a preliminary left multiplication of £ and f by an appro- 
priate matrix will produce the desired standard situation. (See Chapter 25, 
Section 2. ) 

We assxime that E, C, and d are known exactly, or at least that their 
errors are very small relative to those of f. 
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Let C be an X n matrix and let E be X n. We ass^mie that 

mj < n 
mi + m^ > n 
Rank (C) = 

([e]) ' " 

With these assumptions Problem LSE has a unique solution vector and all of 
the solution methods to be discussed apply without the need to consider vinusual 
special cases. 

As a small numerical example to illustrate the computational methods 
to be presented we use the same problem that was used in Chapters 20 - 22. 
(See p. 140). 


[o. 4087 

0. 1593j 

d = 0. 1376 

"0. 4302 

0.3516" 


'0. 6593" 



f = 


.0. 6246 

0.3384. 


.0. 9666. 


The computations described in Chapters 20 - 22 were done using a relative 
-8 

precision of 10 whereas intermediate and final results were rounded to about 
four decimal places for publication. In this note we begin with the published 
intermediate results when applicable and compute using a pocket calculator. 

2. SOLUTION METHOD USING A BASIS OF THE NULL SPACE 

This solution method is described in Chapter 20, pages 134-141. It 
may be s\immarized as follows. 

Apply Householder orthogonal transformations to C from the right to 
reduce C to lower triangular form. Apply these same transformations to E 
from the rignt. Denoting the product of these orthogonal transformations by 
the nxH orthogonal matrix K these operations may be represented by the equation: 
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( 1 ) 


n-m^ 


Solve the following lower triangular system for y^: 

C, y, =d 


Compute: 


f = f - Ej 


^1 


Solve the least squares problem: 


^2 ^2 


= f 


( 2 ) 


Compute: 



To compute the covariance matrix, V, for x^ first compute the covariance 
matrix S for 


S = (E^ Ej) 


(3) 


Then the covariance matrix for 


is 



ro oi 



( 4 ) 


and the covariance matrix for x is 


3 



33-807 


V = KUK*^ (5) 

Consider the numerical example given on pp 140-141. 

In this example 


E 


2 


0. 1714 
.0.0885, 


and 


K = 


-0.9317 
-0. 3632 


-0.3632' 

0.9317 


Thus using Eq. (3) - (5) we obtain 

S = (0.037210)"^ = 26.874 


and 


V 


3.545 -9.094" 
-9.094 23. 33 


( 6 ) 


Note that although Eq. (3) is a valid mathematical definition of S it does not 
represent the most stable way to compute S, If Problem (2) is solved using 
Householder transformations, then one would have an upper triangular matrix 
R such that 


QE 


2 



( 7 ) 


where Q is m^ X m 2 orthogonal. 
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Then, as is described in Chapter 12, one could compute S as 

T 

S = R“^ (R”^) (8) 

3. SOLUTION METHOD USING DIRECT ELIMINATION 

This solution method is described in Chapter 21, pp 144-147. It 
may be summarized as follows. 

Assume column interchanges have been done in the augmented matrix 


C 

E 


if necessary, to assure that the first m^ columns of C are linearly independent. 

Use Gaussian elimination to zero all elements below the diagonal in 
the first m^ columns of 



Solve the least squares problem: 


Solve for in 


E2 x^asf 


Cj Xj = d - C^ X2 


(9) 


( 10 ) 
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Then the solution vector is 


X = 



To compute the covariance matrix of x introduce the m^ x(n - m^) matrix H, 
obtained by solving 


H - 


Then from Eq. (10) we may write 


X 


I 



Hx2 


( 11 ) 


Let S denote the expected value operator. Introduce the mean values 


Xj = ^(Xj) 


and 


^2 ■ ^^^ 2 ^ 

These mean values satisfy Eq. (11), i. e. , 

X, . C[‘ d - HXj 

Subtract Eq. (IZ) from Eq. (11) obtaining 


(12) 


(Xj - X j) = -H(x 2 - X 2 > 


(13) 


from which we may write 


X - X = 






(14) 
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Let W denote the (n - nij) X (n - nij) covariance matrix of x^, which from 
£q. (9) may be defined as 


T • 1 
W = E^) 


( 15 ) 


Then using Eq, (14) the covariance matrix V of x can be written as 



-h' 


HWH^ -HW' 

V = 


w[-H^ l] = 



I 


-WH^ W 


Consider the same numerical example as before, solved by this method. (See 
p. 147). We have 


ro. 1839 



LO. 0949J 


and 


H = = (0.4087)’^(0. 1593) = 0. 38977 


We compute 


W = (eJ = (0.042825)"^ = 23. 351 


and 


V = 


3.548 


L-9. 101 


- 9 . 101 

23. 351. 


Note that Eq. (15) is a valid mathematical definition of W but not a 
recommended computational formula. See the rentark at the end of Sec. 2 for 
suggestions for a more stable way of computing W. 
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4. SOLUTION BY WEIGHTING 

This solution method is described in Chapter 22. It may be summarized 
as follows: 

Suppose the data are scaled so that the elements of largest magnitude 

in the matrices C and E are approximately the sa.ne size. Introduce a scale 

2 

factor, t, such that < is smaller than the working precision. For instance 

-4 -9 

set £ < 10 for Univac single precision arithmetic and £ < 10 for Univac 

double precision. 

Solve the least squares problem 



( 16 ) 


using Householder or Givens orthogonal transformations. 

Solving the problem by either of these methods involves triangulariza- 
tion by left multiplication by an orthogonal matrix O: 


C 


Q 


UE 


d 

£f. 



d 



Hi 

|n - mj 

Hi + "’z ■ " 




1 


Then x is obtained by solving the upper triangular system 



X 



The condition number of this problem is very large (about < however this 
does not affect the accuracy of the solution because of the special structure of 
the matrix and right-side vector. 
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The covariance matrix, V, of x is 

-1 


V ^ 


LcEjJ 


r £ -1 




( 17 ) 


C 

,<El, 


has a large condition number 


Even though the triangular matrix 

its inverse can be computed without numerical difficulty. 

Consider the example used before. The weighted problem to be solved 
(see p. 156) is 


0.4087 
0.4302c 
.0. 6246c 

Since we w'ill be using 4 or 5 place decimal arithmetic we could choose any 
value of C < 10 . The point is that for any two numbers, a and b, of compar- 

able magnitude C should be small enough relative to the computational precision 

2 2 2 

so that the computed value of a + (cb) will just be a . For our numerical 

example we will not assign a specific value to C but will use the computational 

2 2.2 

rule that the computed value of an expression of the form a + (Cb) is a when 
a and b are of the same order of magnitude. 

The data arrays of Eq. (18) can be triangularized by Householder 
transformations to obtain the equivalent problem 


0. 1593 
0. 3516c 
0. 3384C 


0. 1376 ■ 
0. 6593c 
LO. 9666 cJ 


(18) 


'-0.4087 

-0. i5930 ■ 


■-0. 13760' 




"^r 



0 

-0. 20698C 

* 



-0. 80403C 




-*2- 



0 

0 


. 0.43606c 


( 19 ) 
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Solving the nonsinguUr system represented by the first two rows of Eq. (19) 
given the solution vector 


-1. 1774’ 

X = 

. 3. 8846. 

Let R denote the leading 2x2 triangular matrix in £q. (19). We compute 

■-2.4468 1.883U'^' 

R'l = 

. 0 -4.8314C’^ 


tien using Eq. (17) we com|mte the covarian^'c matrix V of x as 

’3.5461 -9.0980" 

-9.0980 23. 342 . 

This computational pirocedure looks peculiar in some ways but it is 
valid. For example the upper left element of R”^, namely -2.4468, is entirely 
lost in the roundoff error when the product R”^ (R*^) i* computed and this 
results in the computed V being singular whereas R’^ was clearly nonsingular. 

This is exactly the right thing to happen, however, since tae covariance 
matrix V for problem LSE should be singular and should not be influenced by 
the upper left element of R"^. 

Close analysis of this weigi.ted method (See Exercise 22.40, p. 157) 
shows that with sufficiently small < this is just a sneaky way of performing the 
direct elimination algorithm treated in Sec. 3 of this note (Chap. 21 of the book), 

5. ONE MORE APPROACH 

Still another way of looking at Problem LSE is presented on pp. 141- 
143. As is noted there we expect that this approach may not have practical 
value but may be of theoretical interest. 
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L*et K be the nxn orthogonal matrix defined in Sec. 2 of this note 
(Chapter 20 of the book). Let K be partitioned as 


Define 


K- [k, Kj ])n 
mj n - mj 


E = (EK 2 ) E 


where the superscript "+'* denotes pseudoinverse. Define 


A 

A = 


A 

LE 


Then, as is proved in Chapter 20, the least squares solution of 


A 

Ax 


Lf. 


(20) 


is the same as the solution of problem LSE: 

Cx = d 
Ex S f 

To compute the covariance matrix of x, regarding x as the solution of 
Eq. (20) , we first write 


X = 
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Assuming the covariance matrix of 



is 



mi m2 


it follows that the covariance matrix, V, of x is 


V 



LO 



( 21 ) 




From Eq. (20. 30) on p. 141 we know that A can be written as 


a'*’ = - K 2 (EK 2 )'^ EC'*' , K 2 (EK 2 )'''] |n 


(22) 




m. 


m- 


Substituting Eq. (22) into Eq. (21) gives 


V 


K2(EK2)^ 



(23) 


From Eq. (1) we have 


EK 


2 


E 


2 


and thus Eq. (23) can be written as 


-s. j-T T 

^ = ^2 ^2 ^2 ^2 

^ T* - 1 T* 

= k^{kJ e^) kJ 

This last expression is identical to the right-side of Eq. (5). Thus we 
obtain the same representation of V as in Sec. 2. 
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